/**
* The computation of 1D advection-diffusion equation
* with explicit scheme u_t + a u_x = b u_{xx}
* author: Zhicheng Yang
*
**/				
# include <stdio.h>
# include <math.h>

# define N 200
# define tol 1.0e-6
# define nu 0.4
# define time 2.0
# define a 1
# define b 0.0475 				/* 这里精确解为u(x,t) = e^{20/19 * x - t} */

void initial(double *u,double dx)
{
	int j;
	FILE *fp;

	for (j=0;j<=N;j++)
		*(u+j)=exp(20.0*j*dx/19);
	fp = fopen("initcase1.txt","w");
	for (j = 0;j < N+1;j++)
		fprintf(fp,"%e\n",*(u+j));
	fclose(fp);
}

void output(double *u)
{
	int j;
	FILE *fp;
	fp=fopen("outputcase1.txt","w");
	for (j=0;j<=N;j++){
		fprintf(fp,"%e\n",*(u+j));
	}
	fclose(fp);
}

/* L^2 and L^{\infty} error */
void outputerror(double *u,double dx,int step)
{
	int    j;
	double err,err2,e;
	double exact[N+1];
	FILE   *fp;

	err = 0.0; err2 = 0.0; e = 0.0;
	for (j=0;j<N+1;j++)
		*(exact+j)=exp(20.0*j*dx/19-time);
	for (j=0;j<N+1;j++){
		e=fabs(*(exact+j)-*(u+j));
		if (e>err)
			err=e;
		err2+=e*e;
	}
	err2=sqrt(err2*dx);

	fp=fopen("exact.txt","w");
	for (j = 0;j<N+1;j++)
		fprintf(fp,"%e\n",*(exact + j));
	fclose(fp);
	
	fp = fopen("case1err.txt","a");
	fprintf(fp,"Case1: u_t + u_x = 0.0475 u_{xx}\n with time = %f, dx=%e, nu= %e.\n\n",time,dx,nu);
	fprintf(fp,"Compute Steps:%d\n",step);
	fprintf(fp,"L^{\\infty} error:%e\n",err);
	fprintf(fp,"L^2 error:%e\n",err2);
	fclose(fp);
}

int main(int argc,int *argv[])
{
	int j,k,step;
	double dt,dx,t,courant;
	double u0[N+1],u[N+1];
	
	dx=1.0/N;dt=nu*dx*dx/b;
	
	t=dt;
	initial(u0,dx);
	step=0;
	
	while (t<time){
		if (t+dt<time)
			t+=dt;
		else{
			dt = time - t;
			t = time;
		}
		courant = a*dt/dx;
		for (j=0;j<=N;j++){
			if (j!=0&&j!=N)
				*(u+j)=*(u0+j) - courant*(*(u0+j+1) - *(u0+j-1))/2 + nu*(*(u0+j-1)-*(u0+j)*2+*(u0+j+1));
			else{
				if (j == 0)
					*(u + j) = exp(-t);
				else
					*(u + j) = exp(20.0/19 - t);
			}
		}
		for (j=0;j<=N;j++)
			*(u0+j)=*(u+j);
		step++;
	}

	printf("\n Step=%d\n",step);
	output(u);
	outputerror(u,dx,step);
	return 0;
} 
